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ABSTRACT 

We investigate the effect of primordial helium on hydrogen reionisation using a hydro- 
dynamical simulation combined with the cosmological radiative transfer code CRASH. 
The radiative transfer simulations are performed in a 35.12ft, -1 comoving Mpc box 
i-^ ■ using a variety of assumptions for the amplitude and power-law extreme-UV (EUV) 

spectral index of the ionising emissivity at z > 6. We use an empirically motivated pre- 
scription for ionising sources which, by design, ensures all of the models are consistent 
with constraints on the Thomson scattering optical depth and the metagalactic hydro- 
gen photo- ionisation rate at z ~ 6. The inclusion of helium slightly delays reionisation 
due to the small number of ionising photons which reionise neutral helium instead of 
hydrogen. However, helium has a significant impact on the thermal state of the IGM 
during hydrogen reionisation. Models with a soft EUV spectral index, a = 3, produce 
IGM temperatures at the mean density at z ~ 6, To ~ 10500 K, which are ~ 20 per 
cent higher compared to models in which helium photo-heating is excluded. Harder 
EUV indices produce even larger IGM temperature boosts by the end of hydrogen 
reionisation. A comparison of these simulations to recent observational estimates of 
the IGM temperature at z ~ 5-6 suggests that hydrogen reionisation was primarily 
driven by population-II stellar sources with a soft EUV index, a ^ 3. We also find 
that faint, as yet undetected galaxies, characterised by a luminosity function with a 
steepening faint-end slope («lf < — 2) and an increasing Lyman continuum escape 
fraction (/ CS c ~ 0.5), are required to reproduce the ionising emissivity used in our sim- 
ulations at z > 6. Finally, we note there is some tension between recent observational 
constraints which indicate the IGM is > 10 per cent neutral by volume z ~ 7, and 
estimates of the ionising emissivity at z = 6 which indicate only 1-3 ionising photons 
are emitted per hydrogen atom over a Hubble time at z = 6. This tension may be 
alleviated by either a lower neutral fraction at z ~ 7 or an IGM which still remains a 
few per cent neutral by volume at z = 6. 

Key words: dark ages, reionisation, first stars - intergalactic medium - cosmology: 
theory - methods: numerical 



1 INTRODUCTION 

The last decade has witnessed the establishment of the two 
key pieces of observational evidence which presently shape 
our empirical understanding of the hydrogen reionisation 
epoch. The first is the Thomson scattering optical depth 
inferred from observations of the cosmic microwave back- 
ground (CMB). This measurement provides a constraint on 
the integrated reionisation history, and is consistent with 
hydrogen reionisation beginning no later than z — 10.6 ± 1.2 
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(Komatsu et al. 2011). The second is the signature of 
Hi Lya absorption in the spectra of high redshift quasars; 
observations of the Gunn & Peterson (1965) trough indi- 
cate that the intergalactic medium (IGM) is largely ionised 
by redshifts less than z ~ 6 (Becker et al. 2001; Fan et al. 
2006). These observational data therefore broadly constrain 
hydrogen reionisation to the redshift range z ~ 6 — 12. 

Despite this progress, however, a detailed determination 
of the timing and extent of hydrogen reionisation, as well as 
the exact nature of the sources responsible for driving this 
process, remains elusive. Due to the integral nature of the 
CMB constraint a wide range of extended reionisation his- 
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tories are compatible with the Thomson scattering optical 
depth measurement. The small neutral hydrogen fractions, 
xm ~ 10~ 4 , at which Lya absorption saturates also leave 
room for alternative interpretations of the quasar data (e.g. 
Songaila 2004; Becker et al. 2007). Furthermore, Mesinger 
(2010) has recently pointed out that even at z ~ 5-6, the 
IGM may still harbour large patches of neutral hydrogen; 
the number of currently known quasar sight-lines is insuf- 
ficient to fully rule out this possibility with intergalactic 
Lya absorption observations alone. 

One route to making further progress is therefore 
developing detailed simulations (e.g. Ciardi et al. 2003; 
Iliev et al. 2007; Trac & Cen 2007; Finlator et al. 2009; 
Aubert & Teyssier 2010; Baek et al. 2010), and semi- 
numerical/analytical models (e.g. Choudhury & Ferrara 
2006; Zahn et al. 2007; Mesinger & Furlanetto 2007; San- 
tos et al. 2010; Shull et al. 2011; Raskutti et al. 2012) which 
can be compared to these data to make inferences about 
the reionisation process. However, most existing numerical 
simulations do not explore the effect of hard, helium ion- 
ising photons on the thermal state of the IGM during hy- 
drogen reionisation (although see e.g. Sokasian et al. 2002; 
Paschos et al. 2007; McQuinn et al. 2009 for treatments 
of He ii reionisation at z ~ 3) . This renders the compari- 
son of these models to measurements of the IGM temper- 
ature at z < 6 problematic. In addition, many numerical 
models significantly over-predict the number of ionising pho- 
tons in the IGM relative to observational constraints on the 
Hi photo-ionisation rate from the Lya forest at z ~ 5-6. 
These data are consistent with ~ 1-3 ionising photons emit- 
ted per hydrogen atom over a Hubble time at z — 6. As 
a result, in order for hydrogen reionisation to complete by 
z — 6 and simultaneously match observational constraints 
from the CMB and the background photo-ionisation rate 
at z < 6, reionisation must be an extended process where 
the ionising emissivity increases at z > 6 (Miralda-Escude 
2003; Meiksin 2005; Bolton & Haehnelt 2007; Haardt & 
Madau 2012; McQuinn et al. 2011). Correctly matching 
these "post-reionisation" constraints therefore has impor- 
tant implications for reionisation and the properties of the 
ionising sources in the early Universe. 

In this work we address these issues using radiative 
transfer simulations of reionisation performed using the code 
CRASH (Ciardi et al. 2001; Maselli et al. 2003, 2009; Parti 
et al. 2011). CRASH is a 3D Monte Carlo based code which 
follows the propagation of ionising photons (from both point 
sources and diffuse radiation) and self-consistently calculates 
the evolution of the gas temperature and ionisation state 
of hydrogen and helium in the IGM. Our approach differs 
from previous studies in two important ways. Firstly, we 
include the effect of helium ionising photons on the progres- 
sion of hydrogen reionisation. This is especially important 
for computing the thermal state of the IGM (e.g. Tittley & 
Meiksin 2007; Cantalupo & Porciani 2011; Pawlik & Schaye 
2011), and it enables us to directly compare our simula- 
tions to recent measurements of the IGM temperature at 
z = 5-6 (Becker et al. 2011; Bolton et al. 2012). Secondly, 
instead of using a numerical sub-grid model for the sources 
of ionising photons, the ionising emissivity in our simula- 
tions is matched to the CMB and Lya forest observational 
constraints by design. The goal of this empirical approach 
is to explore the consequences of satisfying these observa- 



tional constraints for reionisation models from the outset, 
instead of tuning free parameters and/or sub-grid prescrip- 
tions within the simulations. 

The structure of this paper is as follows. We begin in 
Section 2 with a discussion of the empirically motivated 
reionisation models used in our analysis, and continue in 
Section 3 with a description of our numerical simulations. In 
Section 4 we demonstrate that our RT simulations match the 
observational constraints on the Thomson scattering optical 
depth and photo-ionisation rate inferred from the Lya for- 
est at z ~ 6, before going on to discuss in detail the effect 
of including helium on the ionisation and thermal state of 
the IGM in Section 5. We perform a comparison of our sim- 
ulations to the observational data in Section 6 and discuss 
the implications for the the properties of ionising sources 
at z > 6. Finally, we summarise and conclude in Section 
7. An appendix presenting selected numerical convergence 
tests of our simulations is provided at the end of the pa- 
per. Throughout the paper, the following cosmological pa- 
rameters are used: Q.a = 0.74, Q, m — 0.26, fif, = 0.024/i 2 , 
h = 0.72, n s — 0.95 and as = 0.85, where the symbols have 
the usual meaning. 



2 THE REIONISATION HISTORY 

The primary goal of this work is to model the effect of hy- 
drogen and helium ionising photons on the IGM, rather than 
self-consistently modelling star formation and feedback ef- 
fects during reionisation. Rather than use a sub-grid pre- 
scription for modelling the production of ionising photons 
in our simulations, we shall instead adopt an empirically 
motivated approach which satisfies the observational con- 
straints from the CMB and Lya forest at z ~ 6 by design. 
We achieve this by using a simple semi-analytical model to 
initially guide the choice of ionising emissivity within our 
radiative transfer simulations. 

We first define the total comoving hydrogen ionising 
emissivity to be em [s _1 Mpc~ 3 ], where the volume filling 
factor of H n is obtained by solving (e.g. Madau et al. 1999) 



dQmi 
dt 



em 
<n H ) 



■ QhiiChii ^ e \ Hn Qhii(T'). 



(1) 



Here qhii(T) is the case-A recombination coefficient, (n H ) 
is the mean comoving hydrogen number density, (n c ) is the 
mean comoving electron number density, a = (1 + z)^ 1 and 
Chii = {"■Hii)/( n Hii) 2 is the clumping factor of hydrogen 
within the ionised IGM. 

The He in filling factor is modelled in a similar fash- 
ion; the much higher energy photons (> 54.4 eV) required 
to reionise He n mean that this quantity can be decoupled 1 



1 We have, however, ignored the effect of neutral helium on the 
evolution of the Hn filling factor, but the lower number density 
of helium, combined with the higher energy of the He I ionisation 
threshold, mean it will have only a small effect on H I ionisation 
(e.g. Section 5.1). For soft, stellar-like ionising spectra, Hn and 
He II ionisation fronts will furthermore closely trace each other 
during reionisation (Friedrich et al. 2012). Lastly, note that 
He I ionisation is included in our radiative transfer simulations; 
the calculation here guides the choice of ionising emissivity in our 
simulations only. 
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Figure 1. The evolution of the filling factor calculated for four 
of the reionisation models listed in Table 1. The black, red, blue 
and green curves correspond to model £1.2-al.8, £T.2-ctl-3, £ 1.6- 
«3 and £1.2-a3, while the solid and dotted curves display the 
H II and He III filling factors, respectively. 



from Hi reionisation (e.g. Madau et al. 1999). Denning the 
comoving Hen ionising emissivity as CHeii, we then have 



dQtlcIII EHell 



dt 



QhcIiiChcIii - c ^ eln QHcIIl(r), (2) 



where (n Hc > = Y(l - y) _1 (n H )/4, Y = 0.258 is the 
cosmic fraction of helium by mass, (n c )Hcin = ( n ri) + 
2{nHe), ChoIii = {^Hciii)/{^Hcin) 2 and (n c )nu = (nn) + 
2(nHe)QHcin/QHii- For a power-law spectrum with spectral 
index a, eHeii = 4 — "chi. We shall assume T = 2 x 10 4 K 
and adopt time independent clumping factors Chii = 3 and 
ChcIii = 3 in Eqs (1) and (2). Note, however, the assumed 
clumping factor and temperature are used as a guide only, 
and will be computed self-consistently within our radiative 
transfer simulations. 

We next define the redshift evolution of the total co- 
moving hydrogen ionising emissivity as: 



(Songaila & Cowie 2010), while Eq. (4) (Springel & Hern- 
quist 2003) provides a simple parameterisation for the ris- 
ing emissivity at z > 6 (peaking at z = 9) required by the 
Lya forest data (e.g. Bolton & Haehnelt 2007; Pritchard 
et al. 2010). 

We shall consider three different models for the spec- 
tral shape of the ionising emission, all of which achieve 
Hi reionisation by z > 6 (i.e. Qhii ~ 1)- Our reference mod- 
els (£1.2-al.8 and £1.6-a3) assume a — 1.8 and a — 3, while 
a third model (£1.2-al-3) assumes 30 (70) per cent of the 
sources have a = 1 (3). A spectral index of a = 1.8 is typi- 
cal of quasars (Telfer et al. 2002), while a = 3 is consistent 
with star forming galaxies with metallicities close to solar, 
i.e. population-II stellar sources (Leitherer et al. 1999). The 
third model assumes that a fraction of the sources instead 
have rather hard spectra, a = 1, typical of hard quasars or 
population-Ill stars (e.g. Bromm et al. 2001b). 

In Figure 1, the evolution of both Qhii (solid curves) 
and QhoIii (dotted curves) is shown for model £1.2-al.8 
(black curves), £1.2-al-3 (red curves) and £1.6-a3 (blue 
curves). The models are normalised to have similar comov- 
ing hydrogen ionising emissivities at each redshift, ensuring 
that any differences in the reionisation histories are largely 
due to the different EUV spectral indices. For example, 
Hen reionisation is completed (Qfmii ~ 1) progressively 
later in models £1.2-al.8 and £1.6-a3, which have softer 
ionising spectra compared to £1.2-q1-3. 

Finally, in addition to these three reference reionisa- 
tion histories, we shall also consider two further models; 
£1.2-al.8-H which excludes the treatment of helium, and 
£1.2-a3 which results in a late Hi reionisation at z ~ 6. 
We include the latter to explore the possibility that the vol- 
ume weighted neutral fraction in the IGM at z ~ 7 may be 
greater than 10 per cent. Such a substantial neutral fraction 
is suggested by recent observations, which, if confirmed by 
future investigations, may be in tension with models which 
satisfy constraints on the Thomson scattering optical depth 
and the hydrogen photo-ionisation rate (see Section 6 later 
for further details). The parameters for these reionisation 
models are summarised in Table 1. Using these simple emis- 
sivity models, we now turn to describing our cosmological 
radiative transfer simulations. 
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Here £ is a free parameter which sets the amplitude of the 
emissivity, £ = 14/15, £ = 2/3, a is the extreme-UV (EUV) 
power-law spectral index of the sources and Ofb is the spec- 
tral index of the ionising background; we shall assume the 
same value for both. Eq. (3) is consistent with observa- 
tional constraints on the Hi photo-ionisation rate from the 
Lya forest at z < 6 (Bolton & Haehnelt 2007, see also Sec- 
tion 6.2) and the mean free path 2 for Lyman limit photons 



When the mean free path is much smaller than the horizon 



3 NUMERICAL SIMULATIONS 
3.1 Hydrodynamical simulations 

In order to perform our reionisation simulations, we require 
a model for the intergalactic medium. In this work we use 
a hydrodynamical simulation performed in a comoving cu- 
bic box of size 35.12/i^ 1 Mpc. The simulation was performed 
using the parallel smoothed particle hydrodynamics (SPH) 
code GADGET-3, which is an updated version of the publicly 
available code GADGET-2 (Springel 2005) . A total of 2 x 512 3 



scale, £hi °c ThiA^j («b + 3)af — 1 , where Thi is the Hi photo- 
ionisation rate and Ajn is the mean free path of an ionising photon 
at the Lyman limit. For a fixed photo-ionisation rate, a harder 
(softer) EUV spectral index or a smaller (larger) mean free path 
will therefore increase (decrease) the emissivity. 
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Table 1. Summary of the ionising emissivity models used in this 
work. The columns indicate, from left to right, the name of the 
model, the amplitude of the emissivity, £, the assumed EUV spec- 
tral index of the source spectrum, a, and the percentage of sources 
with that spectrum, f a . The final column indicates whether or not 
helium has been included in the simulations. 



Model 


£ 


a 


U [%} 


He 


£1.2-al.8-H 


1.2 


1.8 


100 


No 


£1.2-al.8 


1.2 


1.8 


100 


Yes 


£1.2-ctl-3 


1.2 


1 (3) 


30 (70) 


Yes 


£ 1.6-a3 


1.6 


3 


100 


Yes 


£1.2-q3 


1.2 


3 


100 


Yes 



dark matter and gas particles were followed in the simula- 
tion, yielding a mass per gas particle of 4.15 x 10 6 /i _1 M@. 
Beginning at z = 16, outputs were obtained from the sim- 
ulation at redshift intervals Az = 0.5 until z = 7, and then 
at Az = 0.4 intervals until 2 = 5. Haloes were identified at 
each redshift using a friend-of-friends halo finding algorithm 
with a linking length of 0.2. Star formation was included 
using a simplified prescription which converts all gas par- 
ticles with overdensity A = p/(p) > 10 3 and temperature 
T < 10° K into collisionless stars. Note that because of this 
simple treatment our simulations do not self-consistently 
model star formation and feedback. Instead, as discussed 
in Section 2, we shall model the ionising emissivity during 
reionisation using our empirically motivated prescription. 

The hydrodynamical simulation also includes the photo- 
ionisation and heating of the IGM by a spatially uniform 
ionising background (Haardt & Madau 2001). This model 
assumes the IGM is optically thin, and that the IGM is 
reionised instantaneously at z = 9. Although we shall re- 
compute the IGM ionisation and thermal state with our ra- 
diative transfer (RT) simulations at all redshifts, including 
the UV background in the hydrodynamical simulation at 
z < 9 is nevertheless important for properly modelling the 
gas distribution. The photo-heating significantly reduces the 
clumping factor of the gas in the hydrodynamical simulation 
due to pressure smoothing (Pawlik et al. 2009), and with- 
out this feedback effect the simulation would over-predict 
the gas clumping factor towards the end of reionisation. On 
the other hand, we note that increasing the mass resolution 
of our simulations would increase the clumping factor and 
hence the rate of recombination in the simulations. However, 
we defer a detailed investigation of the clumping factor to a 
future study. It should be noted though that, while the in- 
clusion of a clumping factor assures a better estimate of the 
gas recombination rate, it does not capture all the relevant 
radiative transfer effects, such as self-shielding. 

3.2 Radiative transfer simulations 

Once the hydrodynamical simulation outputs were obtained, 
the gas number densities, n, temperatures, T (at z > 9 only, 
see Section 3.1) and the halo masses, M, were transferred to 
a 128 3 grid for the RT calculations, which are performed as 
a post-process. The gridded densities and temperatures are 
obtained by assigning the particle data to a regular grid us- 
ing the SPH kernel (e.g. Monaghan 1992). The correspond- 
ing grid for the halo masses is obtained by using the cloud- 



in-cell algorithm (Hockney & Eastwood 1988) to assign the 
haloes identified by the friends-of-friends algorithm to a reg- 
ular grid with the same dimensions. 

The RT is followed using the code CRASH (Ciardi et al. 
2001; Maselli et al. 2003, 2009; Parti et al. 2011), which self- 
consistently calculates the evolution of the hydrogen and 
helium ionisation state and the gas temperature. CRASH is 
a Monte Carlo based ray tracing scheme, where the ionis- 
ing radiation and its time varying distribution in space is 
represented by mult i- frequency photon packets which travel 
through the simulation volume. For further details regard- 
ing the radiative transfer implementation we refer the reader 
to the original CRASH papers. For each output i of the hy- 
drodynamical simulation, the RT is followed for a time 
Ut,i = ta(zi+i) — tii(zi), where tn{zi) is the Hubble time 
corresponding to z\ which is the redshift of output i. The 
gas number density is updated at each hydrodynamical sim- 
ulation snapshot, and between two snapshots it is evolved 
as n(x c ,y c ,z c )(z) = n(x c , y c , z c ){zi){\ + z) 3 /(l + z ; ) 3 , where 
(x c ,y c ,z c ) are the coordinates of cell c and z\ > z > Zi+i. 
Although the current implementation of CRASH is able to 
model diffuse radiation without approximations, in this work 
we choose to use the on-the-spot approximation. The infi- 
nite velocity of light approximation is made and a photon 
packet is considered as lost once it has exited the simulation 
box, i.e. we do not use periodic boundary conditions. 

The emission properties of the sources are derived as 
follows. Guided by our semi-analytical calculations in Sec- 
tion 3.1, we assume that the total comoving hydrogen ion- 
ising emissivity at each redshift is given by Eqs. (3) and 
(4). Thus, the total rate of ionising photons emitted at 
each output of the hydrodynamical simulation is given by 
Ni — £Hi(zi)K:om, where Kom is the comoving volume of the 
simulation. The emissivity, Ni, is then distributed among the 
sources according to their gas mass, i.e. = NiMj/M to t,i, 
where j refers to the source and M to t,i is the total gas mass 
of sources at output i. This method of assigning the emissiv- 
ity avoids assuming an escape fraction of ionising photons 
and a star formation efficiency, which are very uncertain pa- 
rameters. Furthermore, as already discussed this empirical 
approach is designed to be consistent with the existing ob- 
servational constraints on the photo-ionisation rate at z ~ 6. 
Depending on the redshift and number of sources, we emit 
10 5 — 10 6 photon packets per source at each i r t,i, correspond- 
ing to a total of ~ 5 x 10 7 - 10 10 photon packets. At z < 8.5 
the total number is always > 10 9 , assuring convergence of 
the results to less than one percent (in relative terms) in the 
ionisation and neutral fraction for all the species, as well as 
the gas temperature (see the appendix for further details). 

The ionisation fraction in the RT simulations is ini- 
tialised to its equilibrium value at 2i n , while the initial gas 
temperatures correspond to those predicted by the hydro- 
dynamical simulation, and remain so until either a cell is 
crossed by a photon packet or at redshifts z < 9. In the lat- 
ter instance, the temperature is held fixed at the z = 9 value, 
prior to the onset of photo-heating in the hydrodynamical 
simulation. Once a cell is crossed by a photon packet, the 
ionisation fraction and gas temperature are then updated 
self-consistently within the radiative transfer calculation. 

We have performed five RT simulations in total in this 
study, using the models summarised in Table 1. In order to 
assess the effect of including helium on the evolution of hy- 
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Figure 2. The Thomson scattering optical depth computed from 
each of our five radiative transfer simulations: Sl.2-al.8-H (cyan 
long dashed line), £1.2~al.8 (black solid). £T.2-al-3 (red dotted), 
£1.6-a3 (blue dashed) and £1.2-a3 (green dotted dashed). The 
shaded area corresponds to the 7-yr WMAP value of 0.088 ± 0.015 
(Komatsu et al. 2011). 



drogen reionisation, in model ifl.2-aT.8-H we include only 
hydrogen with a fraction by mass (number) of 0.742 (0.92). 
Furthermore, in model £ 1.2-oT-3, where there are two pop- 
ulations of ionising sources with different power-law spectra, 
the EUV spectral indices are assigned to sources randomly 
(i.e. no correlation with the halo mass is assumed) to repro- 
duce the correct relative proportions. Note also that in all 
five models the power-law ionising spectra extend to a max- 
imum frequency of ~ 200 eV and that the contribution from 
X-rays is not included. Finally, due to the large number of 
sources present in the box, to reduce the computational time 
we adopt the clustering technique described and tested in 
Pierleoni et al. (in prep). This approach significantly speeds 
up our simulations; for reference, the number of sources in 
the 35. 12ft" 1 Mpc box is reduced from 68 (80597) to 34 
(14112) at z = 15 (8). 



4 EMPIRICAL CALIBRATION OF THE 
REIONISATION SIMULATIONS 

Before proceeding to discuss the results of our simulations 
in detail, we first compare them to the two key observables 
we deliberately calibrate to; the electron scattering optical 
depth and the background photo-ionisation rate at z ~ 6 
inferred from the Lya forest. As mentioned in Section 2, 
our choice for the reionisation histories in the simulations is 
such that these key observational constraints should auto- 
matically be satisfied. 



4.1 The Thomson scattering optical depth 

We first consider the observational constraint on the inte- 
grated reionisation history, in the form of the Thomson scat- 



tering optical depth, r c . In Figure 2 the evolution of r e is 
shown for all five of our RT simulations, together with the 
value measured by the 7-yr WMAP mission, 0.088 ±0.015 (Ko- 
matsu et al. 2011). The optical depth, r e , is calculated from 
the RT simulations as: 



cctt 



n c (t)dt, 



(•») 



where c is the speed of light, <jt = 6.65 x 10" 
the Thomson scattering cross section, n c = nmi + iHeii + 
2riHeiii is the electron number density in units of cm" 3 and 
m is the number density of species i, with i=Hn , Hen and 
He in. Here n c is evaluated directly from the simulations for 
z > z m i n — 5, which is the redshift at which the radiative 
transfer simulations are stopped. At lower redshift, where 
we do not have simulation outputs, we instead calculate n e 
analytically assuming that: (i) the average density equals 
the cosmological mean density; (%%) hydrogen is completely 
ionised; (Hi) XHcii = 1 (^hoIii = 0) for 3 < z < z m i n and 
iEHell = (xHelll = 1) for z < 3. 

The Thomson scattering optical depth calculated in 
this manner has a value of 0.073, 0.095, 0.094, 0.090, 0.081 
for the simulations fl.2-al.8-H, £ 1.2-al.8, £1.2-al-3, £1.6- 
a3 and £ 1.2-a3, respectively. As expected, these values are 
consistent with those measured by the WMAP satellite (Ko- 
matsu et al. 2011). Note, however, that for model £1.2- 
al.8-H we consider only the contribution from hydrogen. 
The inclusion of helium in these models is clearly impor- 
tant, adding an additional r e ~ 0.02 to the total opti- 
cal depth for £1.2-al.8. This is largely because of the ex- 
tra electrons liberated by the reionisation of helium, but 
will also be partly due to the higher IGM temperatures 
which arise from Hen photo- heating; the temperature de- 
pendence of the H n recombination rate, ami oc T~ 0,7 , 
means higher temperatures will produce a slight increase 
in the Hn fraction and hence the electron number density. 



4.2 The background photo-ionisation rate 

The photo-ionisation rates are compared to the observa- 
tional data in Figure 3. This comparison, however, is less 
straightforward for two reasons. Firstly, the photo-ionisation 
rate is not a direct output from our RT simulations, and so 
we must estimate it indirectly by assuming ionisation equi- 
librium in each cell (x c , y c , z c ), such that: 



Thi = aHii(T) 



n c WHii 



■ 7eHl(r)ne, 



(6) 



where ami and 7 e m are the hydrogen recombination and 
collisional ionisation rate in units of cm 3 s" 1 , respectively. 
All the other quantities have their usual meaning. This will 
be a reasonable approximation for most of the cells in our 
simulation volume after they have been reionised, but will 
break down close to reionisation when non-equilibrium ef- 
fects are important. Secondly, the observational constraints 
on the photo-ionisation rate are derived from the Lya ab- 
sorption observed in z ~ 6 quasar spectra (e.g. Fan et al. 
2006; Bolton & Haehnelt 2007; Calverley et al. 2011). The 
transmitted Lya flux at these redshifts preferentially sam- 
ples highly ionised, underdense regions in the IGM, and so 
we must take care to use similar criteria when comparing to 
volume averaged values in the simulations. 
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In the upper panel of Figure 3 the evolution of the vol- 
ume averaged Hi photo-ionisation rate, Fhi, is shown for 
model £1.2-a3. The different curves display Thi for a va- 
riety of different sub-samples drawn from the simulation 
volume. The black solid curve shows the photo-ionisation 
rate for all cells, whereas the dotted red curve displays the 
data for underdense cells (A < 1) only. The remaining three 
curves again show the photo-ionisation rate in underdense 
cells, but now with the additional condition that xm < 1CF 2 
(blue dashed), 1CT 3 ( green dot-dashed) and 10 4 (cyan long 
dashed). These cuts correspond to ~ 0.13, 0.13, 0.84 per 
cent of the total number of cells in the simulation volume 
at z = 14. At z = 6 the percentages are instead 63, 62, 
and 18, respectively. When all cells are included, the evo- 
lution of Fhi rises to a peak at z ~ 8 (following the rising 
emissivity at z > 6 in Eq. 4) but declines toward higher 
redshift. This is because a larger number of neutral cells 
are present toward higher redshifts, lowering the volume av- 
eraged photo-ionisation rate. The average photo-ionisation 
rate is slightly lower if only underdense cells are included 
because the overdense (and hence first to reionise) regions 
are discarded. In other words, the photo-ionisation rates are 
higher in the overdense cells since the ionising radiation is 
correlated with the underlying density field (see also Iliev 
et al. 2008; Mesinger & Furlanetto 2009). 

At z = 6, by which time all the underdense regions 
in the simulation have been reionised, all curves converge 
to a similar value. Note, however, that in the cases where 
cuts in the neutral fraction are also applied, at z > 6 the 
photo-ionisation rate is always higher compared to the aver- 
age for all the underdense cells (red dotted curve). This is in 
part because the averages are, by definition, only for highly 
ionised cells which are assumed to be in ionisation equilib- 
rium. The difference is more pronounced at z > 8, however, 
when the ionised regions probed are the increasingly rare 
ionised bubbles around sources. We thus also expect higher 
photo-ionisation rates because the selected cells are closer 
to the ionising sources. However, these regions are rare and 
so only provide a small contribution to the overall volume 
averaged ionisation rate. 

In the lower panel of Figure 3 the evolution of the vol- 
ume averaged Thi is shown for all five simulations in un- 
derdense cells which are highly ionised only (xm < 10 -4 ). 
Note that this cut most closely represents the regions of the 
IGM from which the photo-ionisation rates are measured at 
z ~ 6 (Bolton & Haehnelt 2007). The redshift evolution of 
Thi is, as might be expected, similar for all models. Model 
£1.2-cv3 typically gives a smaller photo-ionisation rate due 
to the lower normalisation of the emissivity. On the other 
hand, model fl.2-al.8-H always has a slightly lower value of 
Thi compared to the case including helium, £1.2-al.8. Note, 
however, the photoionisation rates are inferred from Eq. (6) 
rather than directly obtained, and so variations in the gas 
temperature and electron number density in this model will 
be partly responsible for this difference. 

Finally, as required, we find that for all models at z = 6 
the photo-ionisation rates are consistent with the observa- 
tional constraints from the Lya forest (Wyithe & Bolton 
2011) and proximity effect (Calverley et al. 2011), repre- 
sented by triangles and stars with error bars in Figure 3, 
respectively. On the other hand, the photo-ionisation rates 
at z — 5 underpredict the observed values by a factor of 
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Figure 3. Upper panel: The redshift evolution of the volume av- 
eraged photo-ionisation rate, I'm, for model f 1.2-q3. The black 
solid curve shows the photo-ionisation rate for all cells, while the 
dotted red curve displays the data for underdense cells (A < 1) 
only. The remaining three curves show the photo-ionisation rate 
in underdense cells, but now with the additional condition that 
x m < 10 -2 (blue dashed), 10 -3 ( green dot-dashed) and 10 
(cyan long dashed). Lower panel: The redshift evolution of the 
volume averaged Thi for cells with overdensity A < 1 and 
Xm < !0~ 4 only. The curves correspond to the models fl.2~etl.8- 
H (cyan long dashed lines), £1.2-0:1.8 (black solid), £1.2-al-3 
(red dotted), £1.6-a3 (blue dashed) and £1.2-a3 (green dotted- 
dashed). In all panels the triangles and stars display respectively 
the observational constraints from the Lya forest (Wyithe & 
Bolton 2011) and the proximity effect (Calverley et al. 2011). 

2-3, despite the fact we have deliberately used an ionising 
emissivity which agrees with these data when assuming a 
mean free path consistent with recent observational mea- 
surements (e.g. Songaila & Cowie 2010). This discrepancy 
may be understood by recalling that Thi oc em Am, where 
Am is the mean free path at the Lyman limit. Assuming a 
power-law slope for the H i column density distribution of 
P = 1.3, Songaila & Cowie (2010) measure Am =- 84 (49) 
comoving Mpc at z = 5 (6). In comparison, our simulation 
volume is 48.7 comoving Mpc on a side. This sets an effec- 
tive upper limit on the mean free path of ionising photons 
in our simulations which is around half the observed value 
at z = 5. Our small simulation box therefore most likely 
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accounts for this apparent discrepancy, and we caution that 
the ionising emissivity in our simulations is underestimated 
at z < 6 as a result. 



5 THE EVOLUTION OF THE IGM 

IONISATION AND THERMAL STATE 

We have found that our simulations are in reasonable agree- 
ment with both the observed Thomson scattering optical 
depth and background photo-ionisation rate at z = 6, giv- 
ing us confidence that we may now explore the implications 
of these models for the ionisation and thermal state of the 
IGM in further detail. 

5.1 The ionisation fraction 

The volume averaged ionisation fractions predicted by the 
RT simulations are displayed in Figure 4, where the upper, 
middle and lower panels refer, respectively, to the evolution 
of the H ii, He n and He in fractions for the models summa- 
rized in Table 1. Reionisation is largely complete by z = 7 in 
all models (i.e. xhi < 0.05), with the exception of £1.2-a3, 
which has an H i fraction of 0.15 at z — 7. 

Although the aim of this study is not to compare the RT 
simulations with the semi-analytic calculations used to guide 
our choice of ionising emissivity, it is interesting to note that 
the numerical models reproduce the semi-analytic results for 
the H ii evolution remarkably well. However, the agreement 
is to some extent a fortunate coincidence; a different as- 
sumption for the hydrogen clumping factor in Hn regions, 
Chii, or IGM temperature in the semi-analytical model 
would worsen the agreement. The agreement between the 
numerical and semi-analytical evolution of s;hcIii is slightly 
poorer, which is indeed most likely due to slightly different 
values for the clumping factor and/or temperature in the 
two approaches. Nevertheless, the general agreement indi- 
cates that semi-analytical approaches are indeed useful for 
quickly exploring parameter space in reionisation models, 
at least in terms of the volume of the IGM which is ionised. 
This is perhaps not too surprising; both calculations are ef- 
fectively just counting ionising photons and recombinations. 
Indeed, "semi-numerical" schemes which additionally follow 
the topology of reionisation are also in relatively good agree- 
ment with the results of full RT calculations (e.g. Zahn et al. 
2011). 

The long dashed cyan curve in the top panel of Figure 4 
compares the £1.2-evl.8-H model, which excludes helium, to 
the corresponding reference run £1.2-al.8. The abundance 
of Hn in 51.2-al.8-H is slightly higher because all of the 
ionising photons (> 13.6 eV) are used to ionise hydrogen. 
The inclusion of helium in model £1.2-al.8 has a small effect 
on the evolution of the neutral hydrogen fraction, as some of 
the hydrogen ionising photons with energies > 24.6 eV are 
now used to reionise Hei. However, the difference between 
xmi in the £1.2-al.8-H and £1.2-al.8 models is never above 
a few per cent. 

The impact of different spectral energy distributions on 
the ionised fractions can be seen by comparing model £1.2- 
al.8 to models £1.2-al-3 and £1.6-a3. Interestingly, for the 
mixed source model £1.2-al-3, all three ionisation fractions 
(Hn, Hen, Hem) are extremely similar to those of model 




Z 



Figure 4. Upper panel: The evolution of the volume averaged 
H II fraction calculated with the radiative transfer simulations 
for models £T.2-al.8-H (long dashed cyan line), £1.2-al.8 (solid 
black lines). £1.2-ctl-3 (dotted red), £1.6-a3 (dashed blue) and 
£1.2-a3 (dot-dashed green). Note the solid black, dotted red, and 
dashed blue lines are almost indistinguishable. The stars display 
the semi-analytic result for model £1.2-al.8 (see Section 2.2 for 
details). Middle panel: As for the upper panel but for the volume 
averaged Hell fraction. Note that in this case model fl.2-al.8- 
H is not present. Lower panel: As for the upper middle but for 
the volume averaged He III fraction. The stars again refer to the 
semi-analytic result for model £T.2-al.8. 



£1.2-al.8. This is because both the comoving emissivity and 
the number of photons with frequencies above the helium 
ionisation thresholds are very similar in the two models. 
Spectra with power-law indices a = 1, 1.8 and 3 have a 
percentage of ionising photons above the He i (He n) ionisa- 
tion threshold, i.e. above 24.6 eV (54.4 eV), of ~ 52 (19.5), 
34 (7.5) and 17 (1.5) per cent, respectively. In the case of 
the models with the softer ionising spectrum (i.e. a = 3), 
JEHeiii is much lower due to the paucity of higher energy pho- 



8 B. Ciardi et at 




-8 -6 -4 -2 -8 -6 -4 -2 -8 -6 -4 -2 -8 -6 -4 -2 

log(x) 

Figure 5. The percentage of cells in the radiative transfer simulations as a function Hi, Hn, Hell and He III fractions (from left to 
right) at 2=14 (upper row), 9 (middle row) and 7 (lower row). The curves in each panel correspond to a different reionisation history: 
fl.2-al.8-H (long dashed cyan, first two columns only), £1.2-al.8 (solid black), £1.2-al-3 (dotted red), £1.6-a3 (dashed blue) and 
£1.2-a3 (dot-dashed green). 



tons. The softer spectrum is also reflected in the evolution 
of KHellj which is very similar to that of zhii- 

Finally, model £1.2-a3 exhibits very similar behaviour 
to that of £1.6-aS because they have the same spectral 
index, but the ionisation fractions at the same redshift 
are smaller due to the lower amplitude of the comov- 
ing emissivity. Note, however, that both of these models 
have EUV spectral indices which are too soft to complete 
Hen reionisation by z ~ 2.5-3 (e.g. Fig. 1). These mod- 
els are therefore likely inconsistent with the He n Lya forest 
data at z ~ 3 (e.g. Shull et al. 2010; Worseck et al. 2011; 
Syphers et al. 2011) unless the ionising background spectral 
shape hardens at z < 6, perhaps due to the increasing contri- 
bution of quasars to the ionising background. For reference, 
the volume averaged ionisation fractions at z = 14, 9, 7 and 
6 are summarised in Table 2. 

A more quantitative representation of the distributions 
of the various ionised fractions is displayed in Figure 5, 
where from left to right the percentage of cells as a function 
of xm, xb.ii, XHell and ZHein are shown for the five reioni- 
sation models at z — 14 (upper row), 9 (middle row) and 7 
(lower row). At the highest redshifts most of the hydrogen is 
in a neutral state, but as the redshift decreases and reionisa- 
tion proceeds the percentage of ionised cells increases for all 
models. During the final stages of reionisation (represented 
here at z = 7), most of the cells will be fully or almost fully 
(xnxi > 0.9) ionised and, as a consequence, the percent- 
age of cells with a lower ionisation fraction decreases again. 
Model £1.2-al.8-H generally has a slightly higher number of 
highly ionised cells compared to the three reference models. 



Table 2. Summary of the volume averaged ionisation fractions 
and temperature within the RT simulations. The columns indi- 
cate, from left to right, the name of the model, the redshift z, the 
volume averaged ionisation fractions xmi, XHell and XHellli and 
the volume averarged temperature T. 



Model 


2 


zmi 


^Hell 


^Helll 


T [K] 




14 


0.045 






918 


£1.2-oT8-H 


9 
7 


0.695 
0.981 






9760 
11047 




6 


0.998 






10224 




14 


0.038 


0.017 


0.023 


820 


£1.2-q1.8 


9 
7 


0.632 
0.960 


0.410 
0.499 


0.238 
0.464 


10464 
16594 




6 


0.993 


0.472 


0.522 


16998 




14 


0.038 


0.019 


0.021 


804 


£1.2-q1-3 


9 
7 


0.618 
0.953 


0.419 
0.531 


0.211 
0.425 


10190 
16565 




6 


0.993 


0.490 


0.504 


17454 




14 


0.039 


0.035 


0.004 


643 


£1.6-a3 


9 
7 


0.627 
0.957 


0.594 
0.894 


0.032 
0.063 


7674 
11425 




6 


0.993 


0.922 


0.070 


11347 




14 


0.029 


0.026 


0.003 


488 


£1.2-a3 


9 
7 


0.481 
0.852 


0.459 
0.807 


0.023 
0.044 


6020 
10643 




6 


0.938 


0.888 


0.050 


11347 
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This is again because helium is absent in this model; all the 
ionising photons are thus absorbed by hydrogen, enabling 
hydrogen reionisation to proceed slightly more quickly. The 
behaviour of models f 1.6-a3 and f 1.2-al-3 is also rather 
similar to f 1.2-al.8, except f 1.6-a3 (f 1.2-al-3) has slightly 
less (more) cells with very small ionised fractions. This is 
because of the softer (harder) ionising spectra which pro- 
duce proportionally more (less) hydrogen ionising photons. 
As noted previously, the He n and He in ionisation fractions 
for f 1.2-al-3 and f 1.2-al.8 show rather similar behaviour, 
while f 1.6-a3 exhibits much smaller Hem fractions due to 
the presence of fewer hard, helium ionising photons. A situ- 
ation similar to model £ 1.6-a3 applies to f 1.2-a3, with the 
difference that the lower emissivity means reionisation is less 
advanced. 

From this analysis it is clear that including intergalactic 
helium and a treatment of multi-frequency radiative transfer 
has a rather small effect on the ionisation state of hydrogen 
during reionisation. However, the hard ionising photons ca- 
pable of ionising helium will also significantly photo-heat the 
IGM. We therefore now turn to consider the effect on the 
thermal state of the IGM at high redshift. 

5.2 The volume averaged temperature 

The redshift evolution of the volume averaged gas tempera- 
ture in our five reionisation models is displayed in the upper 
panel of Figure 6. This quantity will depend on the volume 
of the IGM already reionised at any given redshift, as well 
as the spectral shape of the sources in the simulation and 
whether or not helium photo-heating is included. The first 
point to note is that at early times (z > 10) model fl.2-al.8- 
H has a volume averaged gas temperature which is ~10 per 
cent higher than the corresponding model with helium, f 1.2- 
al.8. This is due to the slightly larger volume of the IGM 
in which hydrogen is photo-ionised and heated compared to 
the other models. This arises from the fact (as discussed 
earlier) that no hydrogen ionising photons are used to ionise 
neutral helium. Note, however, that by z ~ 10 the inclusion 
of He n photo- ionisation results in a higher average temper- 
ature for f 1.2-al.8 compared to fl.2-al.8-H. In addition, 
in the absence of any additional heating from Hen photo- 
ionisation, the temperature for £ 1.2-al.8-H slightly declines 
at 2 < 9 as the IGM cools. 

The volume averaged temperature evolution does not 
exhibit any substantial difference between models £ 1.2-al.8 
and £ 1.2-al-3, which is expected from the very similar be- 
haviour of the ionisation fractions discussed earlier. On the 
other hand, despite having a similar behaviour for the evo- 
lution of the H n filling factor, the softer ionising spectrum 
used by £ 1.6-a3 produces temperatures 20-30 per cent lower 
than £ 1.2-al.8. This is partly because the volume filling 
factor of He in is smaller in this model, but also because 
the softer spectrum results in less energy (and hence photo- 
heating) per photo-ionisation on average. Lastly, for the case 
of £ 1.2-a3, the volume averaged temperature is ~20-25 per 
cent lower compared to model f 1.6-a3 over most of reionisa- 
tion, but converges to a similar temperature by z = 6. This 
is due to the lower ionising emissivity, and hence smaller fill- 
ing factor of ionised hydrogen, used in model f 1.2-a3 which 
delays the completion of hydrogen reionisation to z ~ 6. 

We can also isolate the effect of the source spec- 
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Figure 6. Redshift evolution of the volume averaged tempera- 
ture. The curves correspond to models fl.2-al.8-H (long dashed 
cyan), fl.2-al.8 (solid black), £1.2-al-3 (dotted red), fl.6-a3 
(dashed blue) and £T.2-a3 (dotted-dashed green) respectively. 
Upper panel: The temperature calculated by averaging over all 
the cells in the simulation volume. Middle panel: The tempera- 
ture calculated by averaging over only those cells with xmi > 0.9. 
Lower panel: The temperature calculated by averaging over only 
on those cells with ^hcIII > 0.9. 



trum from the volume filling factor of ionised regions by 
calculating the volume averaged temperature in Hn and 
He in regions only, i.e. in regions with xmi > £min (middle 
panel) and xhcIH > aWn (lower panel), where x m i n — 0.9. 
We have verified that varying our choice of threshold results 
in similar average temperatures as long > 0.1. The 

gas temperature reaches its maximum value in the Hn and 
He in regions at the highest redshift, when only a small per- 
centage of cells (< 1 per cent) in the vicinity of the first 
sources have been reached by ionising photons and there 
has been very little time for the gas to cool. As reionisation 
proceeds, more cells are ionised, but those that have been 
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Figure 7. Contour plots of the distribution of gas temperature against proper number density for models fl.2-ol.8-H, £ 1 .2-o 1 .8, £1.2- 
ol-3 and £1.6-«3 (from left to right). The colour scale corresponds to the percentage of cells within each contour. The rows refer to 
redshift z=14, 9 and 7 (from top to bottom). The dashed vertical lines correspond to the average density in the box. 



ionised earlier start to cool primarily by adiabatic expan- 
sion (for gas close to mean density) and Compton scatter- 
ing. The net result is the average temperature in H n regions 
decreases until z ~ 12 (when ~ 5 per cent of the cells 
have cehti > aWn). At lower redshifts, an increase in the 
number of cells in Hn regions which have also experienced 
He ii photo- heating, combined with the fact that more cells 
are being reionised per unit time with the increasing emis- 
sivity, results in the volume averaged Hn region tempera- 
tures gradually increasing again toward z = 6. Note, how- 
ever, that for model fl.2-al.8-H, where Hen heating is ab- 
sent, the temperature starts to fall again at z < 8 once 
Hi reionisation is complete and the ionising emissivity be- 
gins to decline. 

The behaviour of the volume averaged temperature in 
the He in regions (lower panel) is broadly similar to the case 
for Hn regions, with a high initial temperature followed by 
cooling. However, in this instance the temperature remains 
almost constant at z < 12. Here the effect of cooling is 
offset by the temperature increase due to freshly ionised 
He in regions which continue to grow at z < 6. Finally, note 
that for both the Hn and Hem regions, models f 1.2-al.8 
and f 1.2-al-3 always exhibit higher temperatures compared 
to the other models because of the energy input from hard 
photons during He n photo- heating. This is of particular rel- 



evance when comparisons with observations are made, and 
will be further discussed in Section 6. 



5.3 The IGM temperature-density relation 

The temperatures in the simulations are examined in more 
detail in Figure 7, which displays the distribution of the 
gas temperature versus the proper number density for 
fl.2-al.8-H, fl.2-al.8, f 1.2-al-3 and fl.6-a3 (from left 
to right). From top to bottom, each row displays the 
temperature-density plane at redshift z=14, 9 and 7. For 
reference, the volume averaged temperatures at z=14, 9 and 
7 for all models are given in Table 2. All cases show com- 
mon features. While initially most of the neutral gas lies 
along a cold (~ 25 K) isothermal locus, as reionisation pro- 
ceeds more cells are photo- heated into a second, multi- valued 
grouping at higher temperature. At z = 14, a plume of hot- 
ter gas extending out to T ~ 10 3 K from the cold group- 
ing toward higher densities is clearly apparent; this is due 
to shocked heated gas in the hydrodynamical simulation. 
Towards the end of reionisation, the vast majority of cells 
have reached their maximum temperature, which depends 
primarily on the ionising spectrum adopted. The fact that 
ionisation proceeds at a faster pace in model fl.2-al.8-H 
is reflected by the temperature behaviour: while at z = 7 
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almost all the cells in case fl.2-al.8-H have been reached 
by ionising photons and thus heated up, in the other three 
models many cells are still cold and neutral. 

There is also a significant amount of scatter in the tem- 
perature at fixed density at all redshifts. This scatter arises 
from the different reionisation history of each cell in the 
simulation (i.e. inhomogeneous reionisation) as well as the 
fact that we do not use monochromatic photons, but rather 
a spectral energy distribution which can also be hardened 
by spectral filtering (Abel & Haehnelt 1999). This differs 
significantly from the tight, power- law temperature-density 
relation expected in the optically thin case following reioni- 
sation (Hui & Gnedin 1997). 

There are also some small quantitative differences in 
the slope and amplitude of the temperature-density relation 
T = T A 7 ~\ which are summarised in Table 3. It has been 
noted both observationally (Becker et al. 2007) and theo- 
retically (Bolton et al. 2004; Tittley & Meiksin 2007; Trac 
et al. 2008; Furlanetto & Oh 2009) that the temperature- 
density relation may be multiple valued and inverted fol- 
lowing Hi reionisation. This occurs because voids tend to 
be reionised last and have therefore had less time to cool. 
The theoretical study of Trac et al. (2008) in particular 
found 7 — 1 ~ - 0.2 at the end of reionisation. These au- 
thors used a larger simulation volume (lOOh, -1 Mpc) com- 
pared to this work, but found the strong correlation be- 
tween the density field and redshift of reionisation in these 
models extends down to scales of lh~ Mpc. We find the 
temperature-density relation is indeed very mildly inverted 
(7 - 1 ~ -0.05) for £1.2-al.8-H at z = 14, but it remains 
close to isothermal for all other models at all redshifts. The 
origin of the diffferences between Trac et al. (2008) and 
this work are not clear. One possibility, however, is that 
Trac et al. (2008) used a rather different prescription for the 
source emissivity based on the star formation implementa- 
tion of Trac & Cen (2007). The ionising photon production 
rate in this model is not calibrated to match constraints 
from the Lya forest data, and it therefore rises continuously 
toward lower redshift. This means that the latter stages of 
reionisation occur more rapidly in their simulations com- 
pared to our model. A more rapid end to reionisation could 
potentially explain the more strongly inverted temperature- 
density relation Trac et al. (2008) find; proportionally more 
of the underdense gas will have been reionised and reheated 
close to the end of reionisation. 



6 IMPLICATIONS FOR REIONISATION 
SOURCES 

In this section we now consider the implications our em- 
pirically motivated simulations for reionisation by compar- 
ing them to observational constraints on the IGM temper- 
ature at mean density, the volume averaged neutral hydro- 
gen fraction and recent estimates of the ionising emissivity 
from measurements of the UV galaxy luminosity function at 
4 < z < 8. 



6.1 The thermal state of the IGM at z ~ 5 - 6 

We first compare our simulations to recent measurements of 
the IGM temperature in Figure 8 (see also Raskutti et al. 
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Figure 8. The IGM temperature at mean density, To, at different 
redshifts (see Table 3). The filled symbols refers to the values 
measured by Becker et al. (2011, black square) and Bolton et al. 
(2012, red triangle). The curves refer instead to the simulated 
results from models fl.2-ol.8-H (long dashed cyan), £1.2-al.8 
(solid black), £1.2-al-3 (dotted red), £1.6-a3 (dashed blue) and 
£1.2-a3 (dotted-dashed green). 



2012). Becker et al. (2011) recently presented constraints on 
the thermal state of the IGM based on Lya forest observa- 
tions in the redshift range 2.0 < z < 4.8. Their temperature 
measurement at z = 4.8 is reported as To = 8930 ± 2020 K 
(2a errors) assuming an isothermal temperature-density re- 
lation (7 = 1). This constraint is shown by the black square 
in Figure 8. At higher redshift, z ~ 6, Bolton et al. (2012) 
have measured the temperature of the IGM within ~ 5 
proper Mpc of seven quasars using the Doppler widths of 
Lya absorption lines. They report a line-of-sight averaged 
temperature at the mean density of To ~ 16200 K. Note, 
however, this constraint is complicated by the fact that these 
quasars also reionise the Hen in their vicinity due to their 
hard ionising spectra. Bolton et al. (2012) therefore also pro- 
vided an estimate for the temperature after subtracting the 
expected heating from the local reionisation of He 11 by the 
quasars, To ~ 7100 K, assuming a quasar EUV spectral in- 
dex of a = 1.5. This latter estimate is displayed in Figure 8 
as the red triangle with 95 per cent confidence error bars. 
Lastly, note that this constraint is dependent on the uncer- 
tain amount of He 11 heating expected from the quasars; as- 
suming a harder (softer) EUV spectral index for the quasars 
would lower (raise) this temperature constraint by several 
thousand degrees. 

Keeping this in mind, the curves in Figure 8 display 
the temperature at mean density, To, calculated in cells 
with xmi > 0.99 (see Table 3) in models £1.2-al.8-H 
(long dashed cyan), £1.2-al.8 (solid black), £ 1.2-al-3 (dot- 
ted red), £1.6-a3 (dashed blue) and £1.2-a3 (dotted-dashed 
green). We estimate the temperature from the simulations 
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Table 3. The tcmpcraturc-dcnsity relation of the ionised IGM in our simulations. The columns indicate, from left to right, the name 
of the model, the redshift z, and the best fit power-law parameters for the power-law temperature-density relation T = ToA 7-1 . Here 
To and 7 — 1 arc calculated only in cells with xhii > 0.99 (columns 3 and 5) and £Hcm > 0.99 (columns 4 and 6). See text for further 
details. 



1\ A^.A^I 

Model 


z 




To [K] 




7-1 






x ml > 0.99 XHein > 0.99 


xhii > 0.99 


ZHelll > 0.99 




14 


16744 




-0.0404 




£1.2-al.8-H 


9 
7 


10525 
9293 


- 


0.0116 
0.0313 






6 


8648 




0.0419 


- 




14 


19705 


17027 


-0.0153 


0.0101 


£1.2-al.8 


9 
7 


14272 
15823 


11735 
14080 


0.0381 
0.0367 


0.0567 
0.0559 




6 


15927 


15008 


0.0341 


0.0357 




14 


18885 


16743 


-0.0055 


0.0123 


£1.2-al-3 


9 
7 


14145 
15826 


11796 
12624 


0.0385 
0.0370 


0.0641 
0.0678 




6 


16236 


14970 


0.0351 


0.0455 




14 


14285 


13828 


0.0014 


0.0179 


£1.6-a3 


9 


10058 


10753 


0.0434 


0.0438 


7 


9922 


11511 


0.0554 


0.0483 




6 


9725 


13386 


0.0594 


0.0465 




14 


14035 


13779 


0.0045 


0.0169 


£1.2-a3 


9 
7 


10049 
10649 


11005 
11852 


0.0425 
0.0437 


0.0412 
0.0424 




6 


10468 


13038 


0.0453 


0.0544 



in this manner to ensure any neutral gas which has yet to 
be ionised is excluded; the temperature measurements from 
the Lya absorption measurements only probe highly ionised 
hydrogen. The simulations which have a soft (a = 3) EUV 
spectral index (£1.2-a3 and £ 1.6-a3) as well as the model 
which excludes helium (£1.2-al.8-H) are similar or slightly 
greater than (within ~ 0.02 dex of the 95 per cent confidence 
interval) the measurement obtained by Bolton et al. (2012) 
at z ~ 6. In contrast, the two models with harder spectra 
(£1.2-al.8 and £1.2-al-3) exhibit significantly higher tem- 
peratures due to additional Hen photo-heating. Similarly, 
the Becker et al. (2011) temperature measurement at z = 4.8 
is also much lower than the predicted simulation tempera- 
tures at z — 5 for the harder ionising spectra. Note that the 
agreement would be even worse if the heating contribution 
from X-rays were included in the simulations. 

These results are thus consistent with a predominance 
of sources with relatively soft (a > 3) ionising spectra 
during hydrogen reionisation, and also with an epoch of 
He 11 reionisation (most likely driven by quasars) which was 
not fully underway until lower redshift (e.g. McQuinn et al. 
2009). We therefore conclude that if a population of sources 
with rather hard spectra, such as mini-quasars (Madau et al. 
2004) or population-Ill stars (Bromm et al. 2001b) were re- 
sponsible for reionising hydrogen, their contribution must 
be either (i) sub-dominant at all redshifts or (ii) confined 
predominantly at early times (z > 9), such that there has 
been sufficient time for the IGM temperature to cool and 
doubly ionised helium to recombine by z ~ 6. This is not 
surprising as population-Ill stars are believed to be present 



at z < 9, but, compared to population-II stars, in negligible 
numbers (see e.g. Tornatore et al. 2007; Maio et al. 2010). 
Becker et al. (2012) have also recently pointed out that rel- 
ative metal abundances in the IGM suggest population-II 
stars produced the bulk of hydrogen ionising photons dur- 
ing reionisation. Similarly, although mini-quasars have been 
investigated by a number of authors as possible sources of 
ionising photons, the general agreement is that their contri- 
bution is not dominant (see e.g. Madau et al. 2004; Miralda- 
Escude et al. 2000). In addition, a model in which reionisa- 
tion were dominated by mini-quasars would most likely over- 
predict also the observed soft X-ray background Salvaterra 
et al. (2005). 



6.2 Ionising photon production 

We next compare the ionising emissivity used in our simu- 
lations to observational estimates based on recent measure- 
ments of the galaxy luminosity function at 4 < 2 < 8. For 
this purpose, we compute the ionising emissivity from galax- 
ies using the recent fit to the redshift evolution of the galaxy 
luminosity function presented by Bouwens et al. (2011). We 
assume a spectral energy distribution e„ oc v° for 912A < 
A < 3000A and e v oc v~ s (i.e. a = 3) for A < 912A, with an 
additional factor of six break at the Lyman limit (e.g. Lei- 
therer et al. 1999; Madau et al. 1999). In addition, we adopt 
two different redshift evolutions for the faint-end slope: the 
Bouwens et al. (2011) best fit a L F = -1.84 - 0.05(z - 6), 
and a steeper faint end slope of o?lf = —1.9 — 0.1(2 — 6). 
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These choices are intented to represent the considerable ob- 
servational uncertainty in the faint-end slope. 

The resulting emissivities are displayed as the hatched 
regions in Figure 9, with the results from the two different 
faint end slope evolutions shown in each panel. The cyan and 
orange hatching assume ionising photon escape fractions of 
/esc = 0.2 and / osc = 0.5, while the lower and upper limits 
to the hatching correspond to the emissivity obtained by 
integrating the Bouwens et al. (2011) luminosity function 
fit to a lower magnitude limit of Muv = — 18 and Muv = 
— 10, respectively. These limits roughly correspond to the 
magnitude limit of the observational data and the expected 
magnitude of a galaxy in a halo with virial temperature 2 x 
10 4 K (Trenti et al. 2010), respectively. These are compared 
to the emissivities used in models £1.2-a3 (solid curve) and 
£ 1.6-cv3 (dashed curve). Observational constraints on the 
emissivity at z < 6 (red circles with error bars) derived 
from measurements of the photo-ionisation rate from the 
Lycv forest (Wyithe & Bolton 2011) and mean free path 
(Songaila k, Cowie 2010) are displayed as red circles with 
error bars. Note again, that the models are by construction 
chosen to match these constraints closely. 

In order to match the emissivity in model £1.2-a3 up 
to z = 8, an extrapolation of the faint end of the luminosity 
function to Muv = —10, a high escape fraction / csc = 0.5 
and a slightly steeper faint-end slope than the best fit of 
Bouwens et al. (2011) are required. Faint (and currently un- 
detected) galaxies are thus required to reproduce the ionis- 
ing emissivity in our simulations. Recent theoretical studies 
indicate the faint end slope may indeed steepen at z > 6 
(Trenti et al. 2010; Jaacks et al. 2012). A rather high Ly- 
man continnum escape fraction is also required from these 
faint galaxies. Although impossible to measure directly at 
z > 6, recent observations indicate the escape fraction at 
z ~ 3 is larger than at later times (e.g. Siana et al. 2010). In 
addition, Rauch et al. (2011) have recently presented obser- 
vations of a morphologically disturbed, faint Lya emitting 
galaxy at z = 3.44 which are consistent with a Lyman con- 
tinuum escape fraction of 50 per cent. These authors note 
that such faint, interacting galaxies may be more common 
at higher redshift, where the increasing importance of grav- 
itational interactions and mergers could provide a plausible 
mechanism for such high escape fractions. 

Finally, the emissivity evolution in our simulations is 
such that a halo with a baryon mass Mb = 10 s Mq at z — 14 
produces ~ 5 x 10 52 phot s -1 and ~ 10 50 phot s -1 at z = 6. 
For comparison, the number of ionising photons emitted by 
a halo with baryon mass Mb = Mtot(f2b/fi m ) can be written 
as (see Iliev et al. 2006): 



N 



A/cscjVphotMb 

m p At 
5 x 10 52 phot s~ 

KomJ V 0.5 7 V: 



(7) 



N, 



phot 



5 x 10 3 / V 108 M a 



Mb 



10 7 yr 
~A~T 



where /* is the fraction of baryons which are converted into 
stars, /csc is the escape fraction of ionising photons, iVphot is 
the number of ionising photons per stellar baryon, m p is the 
proton mass and At is the time between two snapshots of 
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Figure 9. The ionising emissivity in models £1.2 — ct3 (solid 
curve) and £1.6 — a3 (dashed curve) compared to observational 
constraints based on Lya forest data at z < 6 (red circles) and es- 
timates of the emissivity from recent constraints on the luminosity 
function of high redshift Lyman break galaxies by Bouwens et al. 
(2011) (hatched regions). Upper panel: Comparison to the emis- 
sivity derived from the best fit redshift evolution of the luminosity 
function at 4 < z < 8 presented by Bouwens et al. (2011) (see text 
for details) with a faint end slope crp = —1.84 — 0.05(z — 6). The 
cyan and orange hatching assume escape fractions of / eS c = 0.2 
and /csc = 0.5, respectively, while the range of the hatched regions 
corresponds to the emissivity obtained by integrating the lumi- 
nosity function to a lower magnitude limit of Mi im = — 10 and 
A^lim = — 18 (upper and lower limit to hatching, respectively). 
Lower panel: As for upper panel, but now assuming a steeper faint 
end slope for the luminosity function, «lf = —1-9 — 0.1(z — 6). 



the hydrodynamical simulation 3 . Typically, A p hot = 5 x 10 3 
and 1 x 10 4 for population-II stars with a Salpeter IMF 
and a top-heavy IMF, respectively (e.g. Iliev et al. 2006). 
The requirement for a large escape fraction (/ OS c ~ 0.5) 
may be therefore relaxed somewhat if the efficiency of ion- 
ising photon production increases toward higher redshift or 
a top-heavy IMF is invoked (see e.g. Bromm et al. 2001a; 
Schneider et al. 2002). However, as noted in the previous 
section, the IGM temperature measurements appear to rule 
out significant reionisation by metal-free stellar populations, 
at least at z < 9. However, as there are a variety of possible 
parameter combinations which could satisfy the emissivity 
required, it is not possible to set a stringent constraint on 
the individual parameters in Eq. (7). 



3 Note that the physically relevant timescale here is actually the 
lifetime of the stellar population. In practice, however, numerical 
simulations assume a uniform emission of ionising photons within 
each At, so that the total number of emitted photons is conserved. 
For a more extensive discussion on Eq. 7 we refer the reader to 
the original paper. 
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6.3 The volume averaged Hi fraction 

Lastly, we compare our simulations to constraints on the 
volume averaged Hi fraction, xhi, in the IGM at z > 6. As 
discussed earlier, the presently available observational data 
remain inconclusive with regard to the redshift evolution of 
xhi- This is largely because almost all the methods used to 
derive xhi are somewhat model dependent and/or are lim- 
ited by the available data. For example, at z — 5.5, studies 
of the transmitted flux in the Lya forest indicate xhi ~ 1CP 4 
(Fan et al. 2006; Becker et al. 2007; Bolton & Haehnelt 2007) 
in the regions where Lya transmission is detected. However, 
Mesinger (2010) has noted that the relatively small number 
of quasar sight-lines which have been analysed, combined 
with the fact that quasars sit in highly biased regions, does 
not preclude an IGM which is still a few per cent neutral by 
volume at z = 5-6; isolated patches of neutral hydrogen may 
still lurk undetected in the diffuse IGM at these redshifts due 
to the inhomogeneous nature of rcionisation (see also Lidz 
et al. 2007). Indeed, taking an (almost) model independent 
approach, McGreer et al. (2011) calculated a conservative 
upper limit of xhi ~ 0.9 from Lya forest data at z ~ 6.1, 
although a subsample of two deep spectra provided a more 
stringent constraint of xhi ~ 0.5. 

Alternative analyses of higher redshift quasar spectra 
also provide variable estimates. An analysis of a putative 
IGM damping wing in a quasar near-zone at z = 6.28 by 
Mesinger & Haiman (2004) yields xhi ^ 0.2. In contrast, 
Maselli et al. (2007) find that the sizes of quasar near- 
zones are consistent with an IGM which is mostly ionized 
at z ~ 6, with xm ^ 0.06. More recently, an analysis of 
the near- zone in the spectrum of the highest redshift quasar 
yet detected was found to be consistent with xhi ~ 0.1 at 
z = 7.085 (Mortlock et al. 2011; Bolton et al. 2011). How- 
ever, all of these observations probe only the neutral frac- 
tion in the vicinity of these quasars, so the interpretation of 
these measurements with respect to the IGM as a whole is 
again hampered by the inhomogeneous nature of reionisa- 
tion (Mesinger & Furlanetto 2008). Lastly, recent measure- 
ments of a rapid decline in the Lya emitter/Lyman break 
galaxy fraction indicate the neutral fraction may be as high 
as xhi ~ 0.5 at z ~ 7 (Schenker et al. 2012; Pentericci 
et al. 2011; Ono et al. 2011). On the other hand, the effect 
of patchy reionisation and galactic outflows on reionisation 
also complicate the use of Lya emitting galaxies as a probe 
of the volume averaged neutral fraction (e.g. Dijkstra et al. 
2011). 

In Figure 10 we present a comparison between the vol- 
ume averaged neutral fraction predicted by our simulations 
and a selection of these measurements. There are two impor- 
tant points to note here. Firstly, all our simulations lie within 
the (admittedly large) region between the lower and upper 
limits at z ~ 6. However, the Mortlock et al. (2011) measure- 
ment appears to exclude all the models with the exception 
of £ 1.2-a3; the neutral fraction in all the other cases is too 
low and as a consequence the emissivity is too high. Recon- 
ciling these models with the Mortlock et al. (2011) neutral 
fraction at z ~ 7.1 would therefore require a lower ionising 
emissivity which then must remain constant or even increase 
weakly toward lower redshift to simultaneously match the 
z = 6 photo-ionisation rate measurements. On the other 
hand, Bolton et al. (2011) note that uncertainties in the 



abundance of high column density systems and the spectral 
shape of the quasar ionising radiation could weaken the up- 
per limit on xhi, so the significance of this difference should 
be treated cautiously. 

The second (related) point is that all four models which 
include helium predict a neutral fraction at z = 6 between 1- 

6 per cent, which lies 1-2 orders of magnitude above the con- 
straints from the Lya forest opacity. This is in stark contrast 
to the conventional interpretation that the IGM is highly 
ionised, xui ~ 10 _ , by z = 6, although this scenario is 
consistent with the conservative estimates of McGreer et al. 
(2011). This result is perhaps not too surprising; numerical 
models which predict a highly ionised IGM at z = 6 typically 
overpredict the photo-ionisation rate or ionising intensity by 
a factor of two or more (e.g. Iliev et al. 2008; Finlator et al. 
2009; Aubert & Teyssier 2010). This implies that when we 
deliberately match the emissivity in our simulations at z = 6 
to the Lya forest data, the IGM is required to have an ap- 
preciable neutral fraction at z ~ 6. A more highly ionised 
IGM by z — 6 may be obtained by adopting an ionising 
emissivity which increases more rapidly than we already as- 
sume at z > 6, but this would still come at the expense of 
not satisfying the z ~ 7 neutral fraction constraint. 

An important caveat, however, is that most reionisation 
models (including this work) do not correctly resolve Ly- 
man limit systems (although see Kohler & Gnedin 2007; Mc- 
Quinn et al. 2011). Lyman limit systems (LLSs) are expected 
to regulate the mean free path of ionising photons once the 
sizes of ionised bubbles exceed the typical separation be- 
tween these optically thick systems (Gnedin & Fan 2006; 
Furlanetto & Mesinger 2009). Since the Hi photo-ionisation 
rate is proportional to the emissivity and the mean free path, 
Thi oc ehiAhi, correctly modelling LLSs is a crucial ingredi- 
ent for simulating the latter stages of reionisation. Although 
our simulations match the observational measurements of 
Thi by design, the mean free path within the simulations is 
not set by LLSs, but rather the remaining patches of neutral 
gas in the IGM which are furthest from the ionising sources 
(in the case of £ 1.2-a3, this is 6 per cent of the IGM by vol- 
ume at z = 6). A mean free path at z — 6 which is instead 
set by LLSs might allow for an emissivity which is consistent 
with the observed constraints on Thi, but at the same time 
have a lower volume averaged neutral fraction due to the 
smaller volume filling factor of these dense optically thick 
systems. 

Note again, however, that the issue of how one could 
then reconcile the large volume averaged neutral fraction 
of xm > 0.1 at z = 7.1 with (i) a low neutral fraction 
of xhi ~ 10~ 4 at z — 6 and (it) an emissivity at z = 6 
equivalent to ~ 1-3 ionising photons emitted per hydrogen 
atom over a Hubble time remains. Since the emissivity must 
increase at z > 6 for reionisation to complete by z = 6 
(Bolton & Haehnelt 2007), either the IGM is more highly 
ionised at z ~ 7 than recent observations suggest, or the 
IGM is still a few per cent neutral by volume at z = 6 
(Mesinger 2010). 

7 SUMMARY AND CONCLUSIONS 

In this work we have investigated the impact of helium 
on hydrogen reionisation using three dimensional, multi- 
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Figure 10. The volume averaged Hi fraction and its evolution 
with redshift. The filled symbols refer to the observational mea- 
surements by Fan et al. (2006, black square), McGreer et al. (2011, 
red triangle) and Mortlock et al. (2011, blue circle, see text for 
further details). The curves display the results of our radiative 
transfer simulations: £1.2-ctl.8-H (long dashed cyan), £T.2-al.8 
(solid black), £1.2-al-3 (dotted red), £1.6-c*3 (dashed blue) and 
£ 1.2-a3 (dotted-dashed green). 

frequency RT simulations. We performed five simulations 
using different models for the amplitude and spectral shape 
of the ionising emissivity during reionisation. By design, all 
our models are consistent with measurements of the Thom- 
son scattering optical depth and the metagalactic hydrogen 
photo-ionisation rate at z ~ 6. This empirical approach en- 
ables us to explore the consequences of satisfying these ob- 
servational constraints for reionisation. The main outcomes 
of this study may be summarised as follows. 

• The evolution of the volume averaged H n fraction, 
shii, is very similar for all models with the same hydro- 
gen ionising emissivity independent of the EUV spectral 
index. However, the spectral energy distribution has a 
strong impact on the evolution on the volume averaged 
Hen and Hem fractions, SHcii and aiHein- Models with a 
soft power-law EUV index, a — 3, produce a much lower 
2'Hcin compared to models in which harder photons are 
present. The inclusion of helium in the RT simulations 
furthermore slightly delays reionisation due to the small 
number of ionising photons which reionise neutral helium 
instead of hydrogen. 

• The choice of EUV spectral index has a significant effect 
on the evolution of the volume averaged IGM temperature 
during reionisation. At z ^1 10, model fl.2-al.8-H (without 
helium) has a volume averaged temperature which is ~ 10 
per cent higher than the corresponding model including 
helium, £1.2-ad.8, due to the slightly larger volume of the 
IGM which is photo-ionised by this time. However, at lower 



redshift the inclusion of He n photo-ionisation results in a 
higher volume averaged temperature for £T.2-al.8. In com- 
parison, despite exhibiting behaviour similar to £T.2-al.8 
and £"1.2-0:1-3 for the evolution of the Hn filling factor, 
the softer ionising spectrum used in £T.6-a3 produces 
volume averaged temperatures which are 20-30 per cent 
lower than £T.2-al.8. This is partly because the volume 
filling factor of Hem is smaller in this model, but also be- 
cause the softer ionising photons produce less photo-heating. 

• The temperature (and ionisation fraction) distributions 
in the simulations exhibit a significant amount of scatter 
at all redshifts. This scatter arises from the different 
reionisation history of each cell in the simulations (i.e. 
inhomogeneous reionisation) as well as the fact that we 
do not use monochromatic photons, but rather a spectral 
energy distribution which can also be hardened by spectral 
filtering. This differs significantly from the tight, power-law 
temperature-density relation expected for an optically 
thin IGM following reionisation. We find the temperature- 
density relation for ionised gas is typically isothermal or 
mildly inverted during hydrogen reionisation. 

• A comparison with recent estimates of the IGM tem- 
perature at z ~ 5 — 6 from Lya absorption in the spectra 
of high redshift quasars suggests that hydrogen reionisation 
is mainly driven by sources with a soft spectral energy 
distribution, a < 3. The simulations with harder spectral 
indices produce temperatures which are larger than the 
observational constraints. We conclude that population-II 
stellar sources are likely to provide most of the ionising 
photons during reionisation, and the spectral shape of the 
ionising background must harden at z < 6 due to the 
increasing importance of quasars if Hen reionisation is to 
complete by z ~ 3. If sources with rather hard spectra, such 
as mini-quasars or population-Ill stars were responsible 
for reionising hydrogen, their contribution must be either 
small or confined to z > 9 to give sufficient time for the 
IGM temperature to cool and for doubly ionised helium to 
recombine by z ~ 6. 

• In order to reproduce the ionising emissivity in our 
simulations at z > 6, we find that the best fit to the 
evolution of the galaxy luminosity function presented by 
Bouwens et al. (2011) at 4 < z < 8 requires extrapolation to 
faint UV magnitudes (Mtjv = —10), as well as a steepening 
faint end slope cvlf < —2 and a high Lyman continuum 
escape fraction / csc = 0.5. Faint, low mass galaxies are 
therefore necessary for providing the required number of 
photons during reionisation, in agreement with several 
other complementary studies. 

• There is some tension between the empirically moti- 
vated ionising emissivity used in our simulations and recent 
observational constraints on the IGM neutral fraction which 
indicate that xm > 0.1 at z ~ 7.1. The ionising emissivity 
inferred from the Lya forest at z = 6 is equivalent to only 1— 
3 ionising photons emitted per hydrogen atom over a Hubble 
time, implying reionisation is extended and that the emis- 
sivity must increase at z > 6 if reionisation is to complete 
by z = 6 (Miralda-Escude 2003; Bolton & Haehnelt 2007). 
However, an increasing emissivity at z > 6 is inconsistent 
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with a large neutral fraction at z ~ 7 in our simulations un- 
less the observations are overestimates or the IGM remains 
a few per cent neutral by volume at z = 6 (see e.g Mesinger 
2010). 

Our results highlight the importance of reproducing post- 
rcionisation constraints such as the IGM temperature and 
background photo-ionisation rate for constraining reionisa- 
tion models. While these simulations were designed mainly 
to investigate the impact of helium on hydrogen reionisation 
and the sources of ionising photons at high redshift, the vol- 
ume used is too small to allow a more detailed discussion 
on helium reionisation (which is thought to be driven by 
quasars and to be complete at z ~ 2.5 — 3) and a more accu- 
rate comparison with observational constraints at z < 6. We 
will postpone this further analysis to a future work, together 
with a more thorough investigation of the impact of unre- 
solved small scale high density peaks. The latter will be par- 
ticularly important for regulating the tail-end of the reioni- 
sation process and for setting the thermal state of the IGM 
by absorbing photons close to the Hi and Hen ionisation 
edges. Including these effects in numerical models is there- 
fore necessary for refining the comparison of simulations 
with observations at z < 6. 
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APPENDIX A: CONVERGENCE TESTS 

As discussed in Section 3.2, depending on the redshift 
and number of sources, we emit 10 s — 10 6 photon pack- 
ets per source at each t r t,i, corresponding to a total of 
~ 5 x 10 7 — 10 10 photon packets. While it is computation- 
ally too expensive to run a full simulation with an order of 
magnitude more photon packets, we have run tests on single 
snapshots and on a limited number of consecutive snapshots 
at high redshift. In Figures Al and A2 the distribution of 
different species and gas temperature, respectively, is shown 
for run £1.2-al.8 (black solid lines) and for the same simu- 
lation with 10 times more photon packets (red dotted). The 
results are shown down to the lowest redshift reached by the 
higher resolution simulation, i.e. z = 10.5, which is obtained 
using 12 snapshots of the hydrodynamic simulation. It is ev- 
ident that an excellent convergence has been reached both 
for the H and He species and the gas temperature, with the 
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Figure Al. The percentage of cells in the radiative transfer simulations as a function of Hi, Hn, Hell and He III fractions (from left 
to right) at z=14 (upper row), 12.5 (middle row) and 10.5 (lower row). The curves in each panel correspond to model £1.2-al.8 (solid 
black lines) and the same model run with 10 times more photon packets (red dotted). 
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exception of cells with xun < 10 -6 and xueii < 10~ 4 . Tests 
using only one snapshot at lower redshifts (i.e. following the 
radiative transfer starting from a non neutral configuration) 
show a similar convergence, but they do not account for 
differences between the two runs which might have accu- 
mulated if the full reionisation history were followed. The 
above Figures though demonstrate that such differences are 
negligible. 
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Figure A2. The percentage of cells in the radiative transfer sim- 
ulations as a function of the gas temperature T at 2=14 (upper 
row), 12.5 (middle row) and 10.5 (lower row). The curves in each 
panel correspond to model £1.2-al.8 (solid black lines) and the 
same model run with 10 times more photon packets (red dotted). 



